FRACALMO PRE-PRINT www.fracalmo.org 
Journal of Computational and Applied Mathematics, 
Vol 207, No 2, pp. 245-257 (2007). 



The role of the Fox- Wright functions 
in fractional sub-diffusion of distributed orderS 

Francesco MAINARDI (1) and Gianni PAGNINI^ 2 ) 

W Department of Physics, University of Bologna, and INFN, 
Via Irnerio 46, 1-40126 Bologna, Italy 
f rancesco . mainardiOunibo . it 

( 2 ) National Agency for New Technologies, Energy and the Environment, 
ENEA, Centre "E. Clementel", 
Via Martiri di Monte Sole 4, 1-40129 Bologna, Italy 
gianni . pagniniObologna. enea. it 

Abstract 

The fundamental solution of the fractional diffusion equation of 
distributed order in time (usually adopted for modelling sub-diffusion 
processes) is obtained based on its Mellin-Barnes integral representation. 
Such solution is proved to be related via a Laplace-type integral to the 
Fox- Wright functions. A series expansion is also provided in order to 
point out the distribution of time-scales related to the distribution of the 
fractional orders. The results of the time fractional diffusion equation of a 
single order are also recalled and then re-obtained from the general theory. 
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1 Introduction 

The Wright function is defined by the series representation, valid in the whole 
complex plane, 

A>-1, ,eC, Z ec. (1.1) 

It is an entire function of order 1/(1 + A), that has been known also as generalised 
Bessel functio^. 

lr This paper is based on an invited talk given by Francesco Mainardi at the international 
conference Special Functions: Asymptotic Analysis and Computation, which took place in 
Santander (Spain) on 4-6 July, 2005. The conference was organized in honor of Nico M. 
Temme, who celebrated his 65-th birthday in 2005. Selected papers presented at the conference 
are published in this special issue of JCAM with the organizers as Guest Editors: Amparo 
Gil, Javier Segura (Universidad dc Cantabria, Santander) and Jose Luis Lopez (Univcrsidad 
Pblica de Navarra, Pamplona). 

2 When A = 1 the Wright function can be expressed in terms of the Bessel function of order 
v = n-l. In fact we have J M _i(z) = (z/2)" -1 Wi,p (~z 2 /A) . 
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Originally, E.M. Wright introduced and investigated this function with the 
restriction A > in a series of notes starting from 1933 in the framework of the 
asymptotic theory of partitions [55J [5(}1 [57J . Only later, in 1940, he considered 
the case — 1 < A < 55]. We note that in the handbook of the Bateman 
Project [T5J (see Vol. 3, Ch. 18), presumably for a misprint, A is restricted to 
be non negative in spite of the fact that the 1940 Wright's paper is cited. 

For the cases A > and — 1 < A < we agree to distinguish the 
corresponding functions by calling them Wright functions of the first and second 
type, respectively. As a matter of fact the two types of functions exhibit a quite 
different asymptotic behavior as it was shown more recently in two relevant 
papers by Wong and Zaho [ 53, 54] . The case A = is trivial since it turns out 
from (1.1) W ^{z) = exp(z)/T(/i). 

Following a former idea of Wright himself [57J , the Wright functions can be 
generalized as follows 



k=0 



(1.2) 



where z 6 C, {a,;, bj} 6 C, {Ai, Bj} £ R with A4, Bj 7^ and i = 1,2, ... ,p, 
j = l,2,...,q. An empty product, when it occurs, is taken to be 1. The 
following alternative notations are commonly used 



(b„B 



(61, -Bi), 



, (a p , A p ) 



(1.3) 



Then, the standard Wright function (1.1), being obtained from (1.2) when p = 
and q = 1 with B% = A > — 1 , 61 = /i, reads 



(PA)' 



(1.4) 



All the above functions are known to belong to the more general class of the 
Fox H functions introduced in 1961 by C.Fox [16]. For more information the 
interested reader is referred to the specialized literature including the books 
[2"5] |5g ] I3T 1 H51 [50]. and the relevant articles g31 [21 123 GO] . In particular, we 
recommend the article by Kilbas et al. [25J where the authors have established 
the conditions for the existence of p ^ q {z), see there §2, and provided its 
representations in terms of Mellin-Barnes integrals, §3, and Fox H functions, 
§4. 

For the sake of reader's convenience we devote the Appendix for a short 
outline of the H functions in order to understand the Fox representation of the 
standard and generalized Wright functions of the first and second type that we 
shall introduce in the following. More appropriately, following (SUEI], we can 
refer to the generalized Wright functions simply to as the Fox- Wright functions. 

The Fox notation for the standard Wright functions depends on their type 
and reads 



W x ^(z) := J2 



k\T(Xk + fi) 



H, 



1.0 

0,2 



H 



1,1 



(0,1)5(1-0, A) 

-;(m>-a) 
(0,1);- 



, A>0: 



1< A < 0. 



(1.5) 
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Putting (bi, Bi) = (fx, A), we have for the generalized Wright function: 



(1.6) 



u p,q+i 


—z 


I H 1,p 
{ u p+i,q 


—z 



(1 



2 j> A?)l>p 



(0,l);(l- M ,A),(l-&,,^) 2i? 



A > ; 



(1.7) 



-1 < A < 0. 



(1 - aj,Aj) 1}P ; 0,-A) 
(0,l);(l-^,S,) 2j9 

In this paper we shall show the key-role of the standard and generalized Wright 
functions of the second type for finding the fundamental solutions of diffusion- 
like equations containing fractional derivatives in time of order /3 < 1. In 
the physical literature, such equations are in general referred to as fractional 
sub- diffusion equations, since they are used as model equations for the kinetic 
description of anomalous diffusion processes of slow type, characterized by a 
sub-linear grow of the variance (the mean squared displacement) with time. 
For an easy introduction to anomalous diffusion and fractional kinetics see the 
popular articles [551 HH] • 

In addition to the simplest case of a single time-fractional derivative, more 
generally we can have a weighted (discrete or continuous) spectrum of time- 
fractional derivatives of distributed order (less than 1): then we speak about 
fractional sub-diffusion of distributed order. We note that only from a few years 
the fractional diffusion equations of distributed order have been investigated, 
over all to describe processes of super-slow diffusion. These processes are 
characterized by a variance growing as a power of the logarithm of time rather 
than as a linear combination of powers with exponent less than 1. 

We shall devote Section 2 to the simplest case of the time- fractional diffusion 
equation of a single order: here we discuss how to obtain the fundamental 
solution that will be expressed in terms of a (standard) Wright function of 
second type. In Section 3 we shall consider the equations of distributed order. 
Starting from a generic distribution of fractional derivatives, we provide some 
representations of the fundamental solution involving Fox- Wright functions 
of the second type. In Section 4, as a check of consistency, we derive the 
fundamental solution for the single order as a particular case of the general 
representations. Finally, the main conclusions are drawn in Section 5. 



2 The time-fractional diffusion equation 
of single order 



It is well known that the fundamental solution (or Green function) of the 
standard diffusion equation 



d d 2 

—u(x,t) = -^u(x,t), ieR, t > 0, 



(2.1) 
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i.e. the solution subjected to the initial condition u(x, 0) = S(x) (the generalized 
Dirac functional) , is the Gaussian probability density function (pdf) 

u(M) = ^*- 1/2 e-* 2 /(4*), (2 . 2) 
that evolves in time with second momenlQ growing linearly with time, 

fi 2 (t) := / x 2 u{x, t) dx = 2t . (2.3) 



We note the scaling property of the Green function, expressed by the equation 

u(x, t) = t~ 1/2 U(x/t 1/2 ) , with U(x) := u{x, 1) . (2.4) 

The function U(x) depending on the single variable x turns out to be an even 
function of x, that is U{x) = U(\x\), and is called the reduced Green function. 
The positive variable X := \x\jt x l 2 is known as the similarity variable. 

By replacing in the standard diffusion equation (2.1) the first-order time 
derivative by an integro- differential operator interpreted as a time fractional 
derivative of order (3 G (0,1], we obtain a generalized diffusion equation, the 
parabolic character of which is preserved. We call it the time-fractional diffusion 
equation of order [3 and, consistently with (2.1), we write it as 

pu(i,i) = pti(i,t), xeR, i>0, 0</3<l, (2.5) 
with u(x, 0) = 5(x). 

In (2.5) dP/dt 13 denotes the fractional derivative (of Caputo type) of order (3, 
whose definition is more easily understood if given in terms of Laplace transform. 
Let f(t) be a sufficiently well-behaved (generalized) function on t > with 
Laplace transform L {f(t); s} = f(s) = J* °° c — s ^ f(t) dt. We have 

L(^f(ty ) s\ = s^f(s)-sf 3 - 1 f(0 + ) with </?<!, (2.6) 



if we define: 



f 1 /"* dftr) dr 

:/(*) for 13 = 1 



dt' 

For < (3 < 1 we can also write the fractional derivative (2.7) in each of the 
following two forms, 



d? 1 



dtP JK ' T(l-/3) dt 



f(r)-f(0 + ) 



dr 



(2.8) 



3 Remar.K: Through this paper we are working formally in that we assume a suitable space 
of generalized functions where it is possible to deal at the same time with delta functions, 
integral transforms of Fourier, Laplace, Mellin type, and fractional integrals and derivatives. 

4 The centred second moment provides the variance usually denoted by cr 2 (t). It is a 
measure for the spatial spread of u(x,t) with time of a random walking particle starting at 
the origin x = 0, pertinent to the solution of the diffusion equation (2.1) with initial condition 
u(x,0) = S(x). The asymptotic behaviour of the variance as t — » oo is relevant to distinguish 
normal diffusion (a 2 (t)/t — * c > 0) from anomalous processes of sub-diffusion (a 2 (t)/t — » 0) 
and of super- diffusion (a 2 (t)/t — > +oo). 
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dy»- /w T(i~f3)dt 



™ dr 



(t-T)/ 3 



We refer to the fractional derivative denned by (2.7) as the Caputo fractional 
derivative, since it was formerly applied by Caputo in the late sixties for 
modelling dissipation effects in Linear Viscoelasticity, see e.g. [SJ O EE [33J ■ The 
reader should observe that Caputo's definition differs from the usual one named 
after Riemann and Liouville, which is given by the first term in the RHS of 
(2.7), see e.g. [H|46]. For more details we refer e.g. to [21 ] l27 l 144] . 

Returning to Eq. (2.5), its fundamental solution can be obtained by applying 
in sequence the Fourier and Laplace transforms to the equation itselQ 

Let f(x) be a sufficiently well-behaved (generalized) function onieR with 
Fourier transform F {f(x); k} — J(k) — f_°° e tKX f[x) dx , n £ R . We have 

F li/W;4 = -« 2 /(«) ( 2 - 10 ) 



dx 2 ' 

and for the Dirac generalized function S(x) we have S(k) = 1 . Then, in the 
Fourier-Laplace domain our Cauchy problem (2.3) appears, after applying the 
formulas (2.6), (2.10), in the form s^ 3 u(k, s) — s' 3-1 = — k 2 u(k, s) , from which 
we obtain 

u( K ,s) = -s o, 0</3<l, 5R(s)>0, k£R. (2.11) 

+ K 

To determine the Green function (that is expected to be symmetric in x) in 
the space-time domain we can follow two alternative strategies related to the 
different order in carrying out the inversion of the Fourier-Laplace transforms 
in (2.11). 

(51) : invert the Fourier transform getting u(x, s) and then invert this Laplace 
transform; 

(52) : invert the Laplace transform getting u(k, t) and then invert this Fourier 
transform. 

Strategy (SI): Recalling the Fourier transform pair, 



?> a „-H6 1/2 



b + K 2 2b 1 / 2 
and setting a = s 13 ^ 1 , b = s@ we get 



6>0, (2.12) 



s /3/2-l | | n/2 

u(x,s) = ^— e-l^l 5 , 0</3<l. (2.13) 



5 The time-fractional diffusion equation was investigated by using Mellin transforms by 
Schneider & Wyss |47l in their pioneering 1989 paper where they adopted the equivalent 
integral form 



u(x, t) = u(x, 0) + 



1 f 

CP) Jo 



d 2 



dr 



{t-rY'l 3 ' 

The time-fractional diffusion equation with the Caputo derivative has been adopted and 
investigated by several authors. From the former contributors let us quote Mainardi, see 
e.g. [301 1311 1321 133) (see also |19l 1201 1221 136] and references therein), who has expressed the 
fundamental solution in terms of a special function (of Wright type) of which he has studied 
the analytical properties and provided plots also for 1 < /3 < 2. 
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The strategy (SI) has been followed by Mainardi [3DJ 13U I3U 133| to obtain the 
Green function in the form 

u(x,t) = t~0 /2 U (\x\/t 0/2 ) , -oo<x<+oo, t>0, (2.14) 

where the variable X :— \x\ft 13 ' 2 acts as similarity variable and the function 
U(x) '■= u(x, 1) denotes the reduced Green function that is expressed in terms 
of a Wright function of the second type. Indeed we have 

U(x) = ~M f (|x|) = ^_ f)1 _ f (-N) , (2.15) 

where the M function of order (3/2 has been introduced and investigated in 
[501 EH E2 133] , see also [Sj . More generally, in the complex plain the function 
Mp (z) is well defined for any (3 € (0, 2) and Vz 6 C by a power series as 

MM) = £ 



fc! r[-/9fc/2+ (1-/3/2)1 
fc r°oo V (2-16) 

= " E HfeT r ^( fc + !)/ 2 ] sin[(7r/3(fc + l)/2] . 

fc=0 

By comparing the power series in (1.1) and (2.16) we recognize that the Mp 
function is indeed a special case of the Wright function of the second type with 
A = — (3/2 and fj, = 1 — (3/2, so that it is an entire function of order 1/(1 — 0/2). 
Noteworthy special cases of this functions are 



M 



(z) = ^= exp (- z 2 /4) , Mi (z) = 3 2/3 Ai (z/3 1/3 ) , (2.17) 



where Ai denotes the Airy function, see e.g. [T| I52j. 

Strategy (S2): Recalling the Laplace transform pair, see e.g. [T51 |2"T| |4"I] . 

& Ep(-dfi) , c>0, (2.18) 

s p + c 

and setting c = n 2 we get 

u(«,t) =E p (-K 2 t (} ), 0</?<l, (2.19) 

where denotes the Mittag-LefHer functional. 

The strategy (S2) has been followed by Gorenflo, Iskenderov & Luchko [T7] 
and by Mainardi, Luchko & Pagnini 35J to obtain the Green functions of the 



6 Let us recall that the Mittag-LefHer function Ep(z) (/3 > 0) is an entire transcendental 
function of order 1//3, defined in the complex plane by the power series 



z 

Eg(z):=y^—, r, 0>O, zSC. 



Originally Mittag-LefHer introduced and investigated (in five notes from f903 to f905) 
this function as an instructive example of entire function that generalises the exponential 
(recovered for (3 = 1). For more details we refer e.g. to [14, 15 18 34 44]. Here we like to 
recall that, for < /3 < 1 and negative argument, Ep preserves the complete monotonicity 
of the exponential: indeed it is represented in terms of a real Laplace transform of a positive 
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more general space-time fractional diffusion equations in terms of Mcllin-Barncs 
integrals. For the time fractional diffusion equation the reduced Green function 
(2.15) now appears in the form: 



U(x) 



i i n +i °° r(i-s) 



7 — lOO 



r(l - p 8 /2) 



x b ds 



cos (kx) Ep (—ft 2 ) dn = 

(2.20) 

with < 7 < 1 and x > 0. We point out that from now on we restrict our 
attention to x > in view of the symmetry of the solution. 

In conclusion, we may represent the solution U(x) given in (2.15) and in 
(2.20) using the general formalism of the Fox- Wright functions, that is in terms 
of a generalized Wright function p ^> q [17], or in terms of a Fox H function [38] . 
as follows 



U{x) = \ * l 



i- 1 - 2 ' 2 > 



2 M 



2 ' 2 / 

(0,1); 



(2.21) 



As proven in [35j we recall that u(x,t) can interpreted as a symmetric spatial 
pd/ evolving in time, with a stretched exponential decay. More precisely, we 
have 



U(x) = ±M § (\x\ 



A/. 



-bx c 



with 



{2tt(2-/3) 2 ' 3 /( 2 -' 3 ) / 3( 2 - 2 ' 3 )/( 2 -«} 

2(2-/3)' v m 

Furthermore the moments (of even order) of u[x, t) are 



A 

2/3-2 



+00 , 
-1/2 



2-/3 



(2.22) 

(2.23) 
(2.24) 



A*2n(i) 



-00 



2n u(x,t)dx= l ( : 2 J l ^ 1 ^ n , n = 0,l,2,...,*>0. (2.25) 
1 (pn + 1) 



Of particular interest is the evolution of the second moment: from (2.25) we 
have 

M*) = 2 * + , 0</3<l, (2.26) 

so that that for /3 < 1 we note a sub-linear growing in time, consistently with 
an anomalous process of slow diffusion (alternatively called sub- diffusion), in 
contrast with the law (2.3) of normal diffusion. Such result can also be obtained 
in a simpler way from the Fourier transform (2.19) noting that 



W»(t) 



8k 2 



u(k = 0, t) . 



(2.27) 



function 

E, 



K Jo 



-at 



t/3-1 



■da , t > , < y9 < 1 . 



cr 2 / 5 +2 0-/3 cos(/3tt) + 1 

but decreases at infinity as a power law with exponent -/3: Ep(-tP) ~ t~' 3 /r(-/3). In 
particular, if j3 = 1/2 we have, for t > and i — > 00, 

-E ; i/2(-v / t) = e* erfc(v / I) ~ l/Cv 7 ^*) , 
where erfc denotes the complementary error function, see e.g. [Tl l52| . 
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3 The time-fractional diffusion equation 
of distributed order 



The fractional diffusion equation (2.5) can be generalized by using the notion of 
fractional derivative of distributed order in tim^|. We now consider the so-called 
time-fractional diffusion equation of distributed order 



b(/3) 



2 r 1 



d0=-^u(x,t), b(0)>O, J b(0)d0 = l, (3.1) 



with x £ R, t > 0. Clearly, some special conditions of regularity and behaviour 
near the boundaries will be required for the weight function b(0). 

Time-fractional diffusion equations of distributed order have recently been 
discussed in [T0 l [TT | [T2 l l48] and in [42]. As usual we consider the initial condition 
u(x,0) = S(x) in order to keep the probability meaning. Indeed, already in the 
paper [10| it was shown that the Green function is non-negative and normalized, 
so allowing interpretation as a density of the probability at time t of a diffusing 
particle to be in the point x. The main interest of those authors was devoted to 
the second moment of the Green function (the displacement variance or mean- 
square displacement) in order to show the sub-diffusive character of the related 
stochastic process by analysing some interesting cases of the weight function 
b{0). 

In this paper we are interested to a more general approach involving a 
generic distribution b((3) in order to provide a general representation of the 
corresponding fundamental solution. By applying in sequence the Fourier and 
Laplace transforms to Eq. (3,1) in analogy with the single-order case, see Eqs. 
(2.5) and (2.11), we obtain, 

= »(*)>0, «GR, (3.2) 

where 

B(s) = { b{0)s^d0. (3.3) 



o 



Before of trying to get the solution in the space-time domain, it is worth to 
outline the expression of its second moment as it can be derived from Eq. (3.2) 
using (2.27). We have 

s) = K 1 " W) + • ■ ■) ' so Ms) = = °' s) = TWA' (3 ' 4) 

Then, from (3.4) we are allowed to derive the asymptotic behaviours of fJ*2(t) 
for t — ► + and t — ► +oo from the asymptotic behaviours of B(s) for s — » oo 
and s — > 0, respectively, in virtue of the Tauberian theorems. 

The expected sub-linear growth with time is shown in the following special 
cases of b{0) treated in [TUl E] . 

The first case is slow diffusion (power-law growth) where 

b(fi) - bxS{fi + b 2 5{(3 - /3 2 ), < fa < 02 < 1, h > 0, b 2 > 0, 6i+6 2 = 1. 



7 We find a former idea of fractional derivative of distributed order in time in the 1969 book 
by Caputo [B] , that was later developed by Caputo himself, see [7] [8] and by Bagley &; Torvik, 
see [2]. 
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In fact 

2 



, - , b 2 T(p 2 + l) 

^ = h l8 ^+i+b 2 ^+i ' S ° M2W ~ \ 2 „ (3 ' 5) 

biTQh + 1 



-t fj2 , t^O, 
I 

-t 01 , t-KX>. 



In [10] . see Eq. (16), the authors were able to provide the analytical expression 
of hi (t) in terms of a 2-parameter Mittag-Leffler function. 

The second case is super-slow diffusion (logarithmic growth) where 

b(/3) = 1, < P < 1 . 

In fact 



logs 



2tlog(l/t), t->0, 



JT 2 (s)=2^— so H 2 (t)~< (3.6) 
S ^- L ) [2 1og(t), t->oo. 

In [TU], see Eqs. (23)-(26), the authors were able to provide the analytical 
expression of (J,2(t) m terms of an exponential integral function. 

Let us now return to Eq. (3.2). Inverting the Laplace transform, in virtue 
of the Titchmarsh theorem we obtain 

«(«,*) = -i f e- rt Im{S(re")} dr , (3.7) 



that requires the expression of — Im {B(s)/{s(B(s) + k 2 ]} along the ray s = r c i7r 
with r > (the branch cut of the functions s@ and s 13 ^ 1 ). By writing 

/ |V N ( p = p(r) = \B(re™)\, 

B [rc l7T ) =pcos(7r7)+ipsin(7r 7 ), < 1 f , ijr ., (3.8) 

V / I 7 = 7M = - arg [B (r e") J , 

after simple calculations we get 

poo —rt 

u(n,t)= K(K,r)dr, (3.9) 



where 

g(«,r) = i ; « 2 PBin(7r7) 

7T K 4 + 2k j /5 COS(7T7) + yCT 

Then, since w(x,i) is symmetric in a;, the inversion formula for the Fourier 
transform yields for x,t > 0, 

i /• +oc f Z" 00 e~ rt 1 

£) = — / cos(nx) < / K(r,K)dr> dn . (3-11) 

n Jo [Jo r J 

To carry out the above Fourier integral we use the method of the Mellin 
transform. Let 

M{f(Z);s} = f*(s)= /(e)e S_1 C 7 i<5R(s)<72, (3.12) 
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be the Mellin transform of a sufficiently well-behaved function /(£) , and let 

i r-f + ioo 

M- i {rw;0 = /(0 = T7 / rwe _a d», (3.13) 



be the inverse Mellin transform, where £ > , 7 = 5R(s) , 71 < 7 < 72 ■ 

Denoting by <-> the juxtaposition of a function /(£) with its Mellin transform 
f*(s) , the Mellin convolution implies 

f°° 1 

h(0 = f(0®9(0--= ~m9(Z/r))dri~h*( S )=r(s)g*(s). (3.14) 
■/o »7 

Then, following [35] (pp. 160-161), we recognize that the Fourier integral in 
(3.11) can be interpreted as a Mellin convolution in k , that is u(x, t) = f(n, t) ® 
g(K,x), if we set (see (3.14) with £ = 1/x, r\ = k) 

roo — rt 

/(«,*) := — — A-(«,r)dr - f*(s,t), (3.15) 

^,x):= J-cosfi) S ) := 9 *(s, X ), (3.16) 

with < 5i(s) < 1. The next step thus consists in computing the Mellin 
transform f*(s, t) of the function /(k, t) and then inverting the product 
f*(s,t) g*(s,x) using (3.16) in the inversion Mellin formula, namely 

1 1 r<?+irx> 

u(x,t) = — — / /*(s,<)T(l - s)sm(Trs/2)x s ds 

1TX Z7TI 7cr— ioo fo 1 7 1 ) 

j 1 ra+ico ^ r(i-s) s w--^; 

= x2ril_ ioc rM r( s /2)r(i- s /2) ^ 

The required Mellin transform /*(s;t) is 

,./ \ / ,00 e- rt fl f 00 K 2 psin(7r 7 ) ,-1,1, ,„ „x 

The term in braces can be computed by making the variable change k 2 — > pfi 
and reads 

s/2+1 1 r sin(^7) . ((s /2 + i)-i 



2p 7T 7 /I 2 + 2/X C0S(7T7) + 1 

(3.19) 

p*/ 2 f r( s /2 + i)r[i-( 8 /2 + i)] - 
2 \ r( 7S /2)r(i- 7S /2) 

where we have used a formula of the Handbook by Marichev, see [39] p. 156, 
Eq. 15 (1), under the condition < 3?(s/2 + 1) < 2, | 7 | < 1. As a consequence 
of (3.18)-(3.19) we finally get 

rC,*) —r — P — l^ + ^-^Vn * , (3.20) 
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Now, using Eqs (3.17) and (3.20) we can finally write the solution as 



1 f°° p-rt 

<x,t) = —J o —F(^ 2 x)dr, 
where F{p 1 / 2 x) is expressed in terms of Mcllin-Barnes integrals: 



(3.21) 



F{p 1 / 2 x) = 



2ni 



(T + 200 



7rr(l - s) 



n , X) r( 7 a/2)r(i - 7 «/2) 

(T+iOO 



{p 1 / 2 x) s ds = 



— / r(l - s)sin(7r 7 s/2)(p 1/2 a;) s ds, 



(3.22) 



and p = p(r), 7 = j(r), we remind it, are related to the distribution b([3) 
according to Eqs (3.3) and (3.8). By solving the Mellin-Barnes integrals by the 
residue theorem we arrive at the series representations in powers of (p 1 / 2 x), 



F(p 1 / 2 x) =np 1 / 2 xJ2 



(-p 1 ' 2 x) k 



,._ fc!r( 7 fc/2 + 7 /2)r(- 7 fc/2 + 1 - 7/2) 
= P 1 ' 2 * E u ^l(k + l)/2) . 



fc=0 



Then, in virtue of Eqs (1.2)-(1.3) we recognize 
F(fi^ 2 x) = np^x o* 2 



(1-7/2,-7/2) (7/2, 7/2)' 



-p x ' 2 x 



(3.24) 



which implies that F{p 1 / 2 x) is a Fox-Wright function of the second type (being 
7 > 0). The Fox representation of the function is 



F(/I 2 x) = nHl; 



1.0 

2 



p x ' 2 x 



— ; (1,7/2) 
(1,1); (1,7/2) 



(3.25) 



In conclusion, the fundamental solution admits the (equivalent) representations: 



u(x, t) 



2 Jo 



3 1/2 0*2 



(l-7/2,-7/2)(7/2,7/2) 



dr, (3.26) 



and 



u(x,t) = - 



00 -rt 

£ o-l,0 



— ; (1,7/2) 
(1,1); (1,7/2) 



dr . 



(3.27) 



If we exchange the order of integration and summation, we have an alternative 
and interesting series representation of the fundamental solution: 



1 ^ {-x) k 



t(x ' i)= 2^ k\ 

k=0 



<Pk{t) , 



(3.28) 



where 

f 

<Pk(t) = / 
Jo 

with p = p(r), 7 = 7 (r). 



■sin[7T7(fc + l)/2] p ( fe+1 )/ 2 dr , 



(3.29) 
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4 The reduction to fractional sub-diffusion 
of a single order v 

In order to check the consistency of the general analysis carried out in the 
previous Section and to explore directions for future work, we find it instructive 
to derive as particular cases the results of Section 2 concerning the fractional 
sub-diffusion of a single order. We now agree to denote this (fixed) order by 
v to be distinguishes from (3 used in the distributed order case (as a variable 
order). This means to consider in Eq (3.1) the particular case 

b{fi) = S{fi-v), 0<i/<1, (4.1) 

so that B{s) = s v and Eq, (3.8) yields 

p = p{r) = r v , 7 = const = v . (4.2) 

In this case the Eqs (3.9)-(3.10) reduce to 

-rt 



f°° c~ rl 
Jo r 

and hence, in virtue of (3.11), 



k p sin(7w) 



7ri \ n Jo « 4 + 2r-, 



2k 2 p cos(7r^) + p 



-rt 



2 



(4.3) 



u(x, t) = — cos(kx) < / — ; — — — dr} dn. 

' / ..i •••"(• u -::r;'i/,-- • | 



With the change of variable r = (jk 2 / u the term in brace reads 

a 2v + 2a v cos(t/7r) + 1 



(4.4) 



^ r * -™ VVt 2 S , , da , (4.5) 



so that 



u(x, t) 



1 f'°° 

= -/ cos{kx)E v {-kH v ) dn, (4.6) 
n Jo 

where E v is the Mittag-Lcfflcr function of order v according to its integral 
representation of footnote ^ . We thus recognize that Eq. (4.6) is consistent 
with Eq. (2.19), and, once applied the scaling relation for the Fourier transform, 
with Eqs. (2.14)-(2.15) and (2.20). 

The consistency with the results expressed in terms of the general formalism 
of Fox- Wright function, that is the comparison between Eqs. (2.21) and (3.26)- 
(3.27), can be obtained in a less direct way because one is required to use the 
scaling relations and the Laplace transform rules of the Fox H functions available 
in the specialized literature. We do not report on this tedious calculations. In 
a more direct and instructive way the consistency with the single-order case 
is shown by using the series representation of the fundamental solution (3.28)- 
(3.29). In this special case the functions tpk{t) turn out to be 

ip k (t) =sinMfc + l)/2] [°° '—l r v{k + l)/2 dr 

Jo r (4j) 
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As a consequence, the solution reads 



u(x, t) = -t-/ 2 ■ - ]T ' mk + l)/2] sm[M(k + l)/2] 

2 n k=o kl (4i 

= l -t~ v ' 2 M, f-^ 
2 2 VW 2 

in agreement with Eqs. (2.14)-(2.16). Of course, only in this special case it is 
possible to single out a common time factor (t~ u / 2 ) from all the functions ipk(t) 
and get a self-similar solution. In general the set of functions tpk(t) give raise 
to a distribution of different time scales related in some way to the distribution 
of the orders of the fractional derivatives. 



5 Conclusions 

The diffusion-like equations containing fractional derivatives in time and/or 
in space are usually adopted to model phenomena of anomalous transport in 
physics, so a detailed study of their solutions is required. Our attention in this 
paper has been focused on the time fractional diffusion equations of distributed 
order less than 1, which are known to be model equations for sub-diffusive 
processes. Specifically, we have worked out how express their fundamental 
solutions in terms of Fox- Wright functions. 

At first we have recalled the main results for the fundamental solution of 
the time fractional diffusion equation of a single order, which arc obtained by 
applying the Fourier-Laplace integral transforms. The required solution turns 
out to be self similar (through a definite space-time scaling relationship), and 
expressed in terms of a special function belonging to the simpler class of the 
Wright functions. Then we were able to adapt the previous techniques for 
obtaining the fundamental solution in the general case of a distributed order. 
For such solution we have provided a representation in terms of a Laplace-type 
integral of a Fox- Wright function, that can be expanded in a series containing 
powers of space and certain functions of time, responsible of the time-scale 
distribution. 

Among the various questions for future research on this topic, particularly 
relevant in our opinion is the possibility to use our analytical results for 
plotting the fundamental solutions in some noteworthy cases of fractional order 
distribution, as it was for the simplest case of a single order. 
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Appendix: The Fox H functions 



According to a standard notation the Fox H function is denned as 

W = i / K^ n ( S )z s ds, (Al) 



where £ is a suitable path in the complex plane C to be disposed later, z s = 
exp{s(log |z| -Margz)}, and 

m n 

A(s) = 11 TQbj - B jS ) , B(s) = l[ T(l a, + A jS ) , (A3) 
i=i i=l 

C(s)= f[ ni-bj + Bjs), D(s)= f[ rfa-Ajs). (A4) 

j—m+l j — n+1 

with < n < p, 1 < m < q, {aj, bj} e C, {Aj, B 3 } £ R+ . An empty 
product, when it occurs, is taken to be one so 

n = B(s) = 1 , m = q<==>C(s) = l, n = p <==> D(s) = 1 . 

Due to the occurrence of the factor z s in the integrand of (A.l), the H function 
is, in general, multi-valued, but it can be made one-valued on the Riemann 
surface of log z by choosing a proper branch. We also note that when the A's 
and B's are equal to 1, we obtain the Meijer's G-functions G™^ n (z). 

The above integral representation of the H functions, by involving products 
and ratios of Gamma functions, is known to be of Mellin-Barnes integral typ^l. 
A compact notation is usually adopted for (Al) : 



Tjm,nf \ rrm,n 

p,q \ ' P,q 



(a,j. Aj)j—i^ n ; (aj, Aj)j— n ^i t p 
(bj, Bj)j = i m ; (bj, Bj)j =m +i.q 



(A5) 



Thus, the singular points of the kernel 7i are the poles of the Gamma functions 
entering the expressions of A(s) and B(s), that we assume do not coincide. 
Denoting by V{A) and V(B) the sets of these poles, we write V(A) C\V{B) = . 
The conditions for the existence of the //-functions can be made by inspecting 
the convergence of the integral (A.l), which can depend on the selection of the 
contour C and on certain relations between the parameters {a^, A} (i = 1, ■ ■ • , p) 
and {bj,Bj} (j — l,...,q). For the analysis of the general case we refer to 
the specialized treatises on H functions, e.g. [101 EH [SD] and, in particular to 
the paper by Braaksma where an exhaustive discussion on the asymptotic 
expansions and analytical continuation of these functions is found, see also |24j . 



8 As historical note we point out that the names refer to the two authors, who in the first 
1910's developed the theory of these integrals using them for a complete integration of the 
hypergeometric differential equation. However, these integrals were first used in 1888 by S. 
Pincherle, see e.g. |37j . Recent treatises on Mellin-Barnes integrals are those by Marichev [39) 
and Paris & Kaminski 1431 . 
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